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l. INTRODUCTION 


A. GENERAL 

In this thesis, several numerical methods were used to analyze the flow 
about one and two airfoils performing unsteady motion in an inviscid, 
incompressible fluid. First, the unsteady motion for a single airfoil was studied 
using the unsteady panel code, U2DIIF, by Teng [Ref. 1]. The primary purpose 
was to verify the phase lag relationships between the airfoil’s motion and the 
build up of aerodynamic forces. To accomplish this, the time dependent output 
of the panel code was converted to harmonic output using a curve-fit algorithm. 
Furthermore, an extensive study on the production of thrust and associated 
efficiency for oscillating airfoils was completed using the U2DIIF code. Finally, 
a theory was presented that related quasi-steady ramp motion to a purely 
steady state phenomenon. 

Next, the vortex interaction between two airfoils was studied using a 
computer code, USPOTF2, developed by Pang [Ref. 12] for unsteady 
incompressible flow. Extensive modifications were made to increase the 
program capabilities. Comparisons were made to flat plate, linear theory using 
the modified code. Again, an analysis of the propulsive forces and efficiencies 
associated with two airfoils was completed using the modified code and a 


similar curve-fit algorithm developed to convert the time dependent output to 


harmonic output. This led up to a systematic study involving the influence of 
an oscillating airfoil in the vicinity of a stationary airfoil. 

Finally, smoke flow visualization experiments were conducted using the 
available facilities located at the Naval Postgraduate School. These 
experiments were undertaken in an attempt to verify the two-dimensional 
theory presented here on the production of thrust associated with an oscillating 
airfoil. Furthermore, a visualization experiment was conducted to better 
understand the influence of an oscillating airfoil in the vicinity of a stationary 


airfoil at high angle-of-attack. 


B. SCOPE 

Chapter II contains the complete single airfoil, theoretical analysis using 
the panel code, U2DIIF. This chapter begins with the development of harmonic 
motion, includes an extensive verification with existing two-dimensional, linear 
theory, and concludes with a presented quasi-steady theory. Chapter Ill 
describes the complete two-airfoil, theoretical analysis. The development of 
the modified version of USPOTF2A is shown along with a complete description 
of the changes. Again, this code is verified with existing two-dimensional, 
linear theory, and used in a systematic study of interacting airfoils. Chapter IV 
and V incorporate some details of the flow visualization experiments and 


selected flow visuzalization photos of the results for the thrust and enhanced 


lift investigation, respectively. Chapter VI contains the conclusion and 


recommendations. 


ll. SINGLE AIRFOIL ANALYSIS 


A. HARMONIC MOTION 


1. Introduction 

The study of two dimensional, unsteady, harmonic motion involves 
an airfoil undergoing sinusoidal pitch or plunge oscillations. As the airfoil 
performs this oscillation, a complex flow field develops as shown in Figure 2.1. 
As the airfoil increases its angle-of-attack, the pressure field around the airfoil 
changes. These changes create a disturbance in the boundary layer due to the 
presence of viscosity. This disturbance builds up on the airfoil surface as it 
travels to the trailing edge in the form of a vortex. The vortex then sheds into 
the medium with a circulation strength equal in magnitude to the increase in 
circulation about the airfoil, but opposite in direction. These disturbances are 
stored in the fluid because the shed vortices convect downstream at the local 
flow field velocity. The counter-rotating vortices induce a sinusoidal flow field 
which further changes the net lift. The result of this complex flow field is a 
time difference or delay in the airfoil’s motion and the induced aerodynamic 


forces. This time difference is known as the phase lag, w [Ref. 1]. 
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Figure 2.1 Oscillating Flow Field 


To simplify the calculations for this type of motion it is common to describe the 
airfoil position and the associated aerodynamic forces with complex variables. 
For pure plunge oscillations the vertical motion of the airfoil is described by the 


real part of the following equation: 


h(t) = he”! (2.1) 
where, hh, is a complex number, and w is the frequency of oscillation. 


Similarly, the lift or moment is described by: 


= res" (2.2) 


where, L, is the theoretical quasi-steady lift given by the expression: 


(2.3) 


This is termed quasi-steady because the angle of attack is represented by h/U. 
The values of r and w represent the magnitude and phase, respectively, of the 
true instantaneous lift relative to the quasi-steady lift. The variables r and win 
general depend on the reduced frequency k, the Mach number M, and the 
Reynolds number. For a nonviscous, incompressible fluid the values of r and 
y will be functions only of k. A complete solution for the oscillating flat-plate 
airfoil in incompressible flow has been obtained by Kussner and Theodorsen and 


is reproduced from reference 3 in Figure 2.2. 
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Figure 2.2. L/Lo diagram 


Harmonic motion is usually characterized by non-dimensional 
variables. The amplitude for plunge motion is made non-dimensional by dividing 
by a reference length, usually the chord or half chord. Pitch amplitude is 
expressed in radians. The nondimensional frequency is termed the "reduced 


frequency” and is given by the following expression: 


k- oF (2.4) 
U 
where, w = frequency of oscillation, b= chord (or half-chord), U= free- 


stream velocity. 


2. Panel Code U2DIIF 

The computer code U2DIIF was developed by Teng [Ref. 2] for the 
study of unsteady, inviscid and incompressible flow over a single airfoil. The 
code was based on a technique called Pane/ Methods developed by Hess and 
Smith [Ref. 4] for steady potential flow problems. This method involves dividing 
the airfoil into many segments or panels. A uniform source and vorticity 
distribution is placed on each panel. In this code, the source strength is 
allowed to vary from panel to panel while the vorticity strength is heid 
constant. This method is based on the fact that the singularity distributions 
automatically satisfy Laplace’s equation, which is the governing equation for 
inviscid, incompressible flow. Furthermore, since the superposition principle 
applies to the Laplace equation, one can build complicated flow fields by a 


combination of simple flows. The unsteady potential flow model is based on 


a method by Basu and Hancock [Ref. 5]. This unsteady model is governed by 
the Helmholtz vortex theorem which requires that any change in circulation 
around the airfoil must be matched by an equal and opposite change in vorticity 
in the wake. This is Known as the vortex shedding process, which is modeled 
in this code by a vortex panel that is shed into the wake at each time step. 
The introduction of a wake model creates a non-linear, unsteady flow problem 
which is solved through an iterative process. No attempt is made here to 
reproduce the work by Teng or to explain the operation of the U2DIIF code, but 
instead the reader is encouraged to review reference 2. One limitation in the 
U2DIIF code which should be noted ts the sensitivity to panel density and airfoil 
thickness. It can be shown that for thin airfoils of less than 8% thickness an 
increasing number of panels (< 50) is required to accurately capture the leading 
edge suction peak. Furthermore, as the airfoil approaches a flat plate there is 
no amount of panels that will capture this peak. Although flat plates or 
NACAOO01 airfoils are used in the following studies, it should be noted that the 
results of this code for airfoils of less than 5% thickness should be considered 
suspect unless it can be validated by other means. 

The code was verified by comparing the time dependent output for 
harmonic motion to experimental results obtained by Giesing [Ref 6.]. In order 
to further verify the U2DIIF code and to present the results for harmonic motion 
in a more useful manner, a program was written to convert the time dependent 


output of lift and moment histories to harmonic output using an iterative, curve- 


fit algorithm (appendix A). Figure 2.3 shows an example of the time dependent 
output, the curve fit solution, and the associated phase angles between the 
motion and the aerodynamic forces. In order to convert the time dependent 
data to harmonic output each run must proceed for at least two cycles. The 
first cycle is discarded because of the transient effects associated with going 
from steady state to unsteady motion as seen in Figure 2.3. The second cycle 
of data (lift and moment histories) is then curve-fitted to the following 


expression: 


F(t) = Amp*Sin(wt + >) (2.5) 


where, Amp = amplitude of motion, w = frequency of motion, and @ = 


phase angle between the motion and the aerodynamic forces. 
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NACAOO12 in pure pitch oscillations at a 
reduced frequency of 2, amplitude of 1°, for 2 cycles. 
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3. Comparison to Flat Plate Theory 
Theoretical results for harmonic motion of a flat plate along with 
some experimental data were published by Halfman and reproduced in 
reference 7. The theoretical results for a flat plate oscillating in two degrees 
of freedom, h and a, were used to compare the U2DIIF code with harmonic 
output. The expressions of aerodynamic force and moment corresponding to 


the harmonic motion given in reference 7 are, respectively: 





L 2 2 _Kwt+,) yf 
— =-,/jA, +i, e@ : ; = tan ‘—— ; 
4qb yy +l o, A, (2.6) 
M 2 2 _ilwttd,) -4 ly 
= /FAy, + ly e “sd yy, =-«s tan" —— (2.7) 
Aqb? M M M Ry 


Furthermore, the coordinate system used in the Halfman experiments was 
defined as positive lift acting down. A direct comparison of the U2DIIF 
harmonic output to these results was made by multiplying the U2DIIF output 
by 2 and adding 180° to the phase angle. This study was conducted with a 
panel density of 100 to accurately capture the suction peak for this very thin 
airfoil. The results are shown in Figures 2.4 through 2.7. It can be seen that 
there is very good agreement between the theoretical results and the U2DIIF 


code. 


11 


Lift Coefficient VS. Reduced FPrequency 

















0.6 
= 
SB oO.5 
= 
ed 
3 
= 
0.4 
0.3 
0.2 
¢ Code 
o Flat Plate Theory 
0.1 Mie os i = 
0.0 =i; rT = | 
0.0 O.1 O.z 0.3 0.4 0.5 0.6 
Reduced Frequency (k) 
Phase Angle VS. Reduced Frequency 
220 = ss iat a <<< = = ] 
=e 200 
o 
bo 
ed 
© 
<o 
—™ 180 
ao 
ra 
160 
140 -_— - 
¢ Code 
a Flat Plate Theory 
120) >) = —— a 
100 Fe uae (ea ei | eS) | pe (oT en ae il 
0.0 o.! 0-2 0.3 0.4 0.5 0.6 


Reduced Frequency (k) 


Figure 2.4. Lift magnitude and phase comparison for a flat plate pitching 
about .37c at an amplitude of 6.7. 
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Moment Coefficient VS. Reduced Frequency 
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| Figure 2.5. Moment magnitude and phase comparison for a flat plate pitching 
about the .37c at an amplitude of 6.7°. 
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Figure 2.6. Lift magnitude and phase comparison for a flat plate in pure 
translation at a magnitude of .0833c. 
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Moment Coefficient VS. Reduced Frequency 
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Figure 2.7. Moment magnitude and phase comparison for a flat plate in pure 
translation at an amplitude of .0833c. 
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B. PROPULSIVE EFFICIENCY 


1. Introduction 

Propulsion or thrust force from harmonically oscillating airfoils is a 
well Known phenomenon. It has long been observed in nature from the low 
speed flapping of a bird’s wing to the high speed flapping created by flying 
insects. The purpose of this study was to obtain the propulsive efficiencies 
associated with pure pitch and plunge harmonic motion using the nonlinear, 
unsteady panel code, U2DIIF. A comparison was made to the results obtained 
by Bosch [Ref. 8] for a flat plate undergoing harmonic motion. Bosch used a 
linear, analytical method to obtain the aerodynamic forces and propulsive 
efficiencies for a flat plate undergoing pure pitch and pure plunge harmonic 
motion. 

The propulsive forces associated with an oscillating foil have been 
experimentally measured in reference 18. In this reference, Scherer measured 
the propulsive forces and moments associated with a rigid foil of finite span 
undergoing large amplitude pitching, and plunging oscillations in water. This 
work was undertaken for the preliminary design of an oscillating foil propulsor 
suitable for use on a shallow-draft boat, such as a ‘ski-barge.’ The results are 
expressed in coefficient form, and may be used in a comparative analysis with 
the theory presented here using the appropriate three-dimensional corrections. 


Sufficient time did not permit a detailed comparison in this paper. 
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2. Theory 

In an inviscid, incompressible flow field propulsive forces develop over 
a harmonically oscillating airfoil. To better understand this process, Figure 2.8 
shows the forces that develop on a harmonically plunging airfoil. It can be seen 
that when the airfoil is moving vertically up or down the motion creates an 
induced velocity component. The force that develops remains perpendicular to 
the relative velocity. This allows for a forward or thrust component of force 
while the vertical component changes sign. This is an over-simplified 
explanation of a complex flow field. The actual flow field that develops will 
depend on the type of motion, the frequency of oscillation, and the intricate 
system of shedding vortices that store the kinetic energy of the motion in the 


fluid. 


Vtotal 


Vtotal 





Figure 2.8. Propulsive force on plunging airfoil. 
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To understand how the mechanical energy of the motion is 
transferred to a propulsive force, it is necessary to develop the relationship for 
the required input work to the airfoil for both pitching and plunging motion. 


First, for a plunging airfoil the average work ts defined as: 


: 
—_ 1 ° 
W = —]| Re[L]*Re[h] at (2.8) 
rl [L] « Re[ A) 
h= het, h = inhe*t, T= 22 ae 
® 


The equation becomes, 


(2.10) 


2n 
W-— f Re[(L_+iL)|[coswt+ isinet] Rel iwh, (coswt+ isinwh | at 
0 


= 


After evaluation the equation reduces to, 


W = Tall (2.11) 


a 
= 
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where, the imaginary lift and reduced frequency are given by, 
L,=pU2bImMC), k= ae (2.12) 


Finally, the expression for average work reduces to, 


W = 7 © Usb kIM(C) (2.13) 


This equation for average work agrees with reference 8 if we set the value 
h,/b =1.0. This equation shows the important functional relationships between 
the input variables of plunge amplitude, reduced frequency, and generated lift 
to the output variable of average work. Furthermore, it shows that the average 
work is proportional to the reduced frequency and the imaginary or ‘out-of- 
phase’ component of lift. 

The development for the relationship of average work for a pitching 


airfoil proceeds in a similar manner starting with the relation, 


; 
W = = [ Re| M] + Re[é) dt (2.14) 
0 


and the final result is, 
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W = ap V2 bk Im(C,) (2.15) 
This result again matches the equation in reference 8 if we set the value of 
Q=1.0 radian. This equation shows the functional relationship between the 
input variables of pitch amplitude, reduced frequency, and the imaginary 
component of moment with the output variable of average work. 
A measure of the efficiency with which work input is transferred into 
propulsive power Is given by the following expression, 


er eA (2.16) 
W 


where, 77 = efficiency, D= average drag or propulsive force, v = velocity. 

This equation for efficiency is simply the drag (propulsive force) times the 
velocity divided by the expression for average work per unit time which is the 
average power output over average power required for the airfoil motion. It is 
clear the efficiency is primarily a function of the thrust that develops for a 
particular type of motion. Since the magnitude of thrust is proportional to the 
induced velocity ahead of the airfoil, the motion which provides the greatest 
induced velocity for the same work input will also provide the greatest thrust 
at the highest efficiency. The two types of motion studied in this paper are 
pitch, and plunge. Pitch motion induces a flow field in front of the airfoil that 


is proportional to the distance the pitch axis is from the airfoil. Furthermore, 


plunge motion can be thought of as pitch motion with the pitch axis at ©. 
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This can be visualized as the motion of a fan blade with pitch motion similar to 
the blade hub, and plunge motion similar to the blade tip. Therefore, it should 
be no surprise that pure plunge motion provides more efficient propulsion than 


pure pitch. 


3. Comparison to Flat Plate Theory. 

Analytical results for a flat plate airfoil undergoing pure pitch and pure 
plunge motion for both propulsive force and efficiency were computed by 
Bosch in reference 8. To compare these results with the nonlinear code, 
U2DIIF, it was necessary to convert the efficiencies to aerodynamic forces. 
Furthermore, Bosch uses a reference length of b/2 (1/2 chord length) where the 
U2DIIF code uses b (chord length) as the reference length. Also, Bosch 


defines the average drag as follows, 


D = ply*C, (2.17) 
where |=1/2 chord. Combining equations 13, 15, and 17 into 16 gives the 


follwing relations: 


Co 2Cp 


= cc (2.18) 
"oten ~ “kim(C,) ’ "purge ~ “kim(C)) 


These equations for efficiency were used to compare the results obtained by 
Bosch. The U2DIIF code was executed at twice the reduced frequency as the 


one desired for comparison due to the different reference lengths. 


21 


In this study there were several problems associated with comparing 
analytical, linear results using a nonlinear code. First, the amplitudes used by 
Bosch: h/b =1.0, and a,=1.0 radian, (used for convenience) were well out of 
the linear range. For this reason much smaller amplitudes were used for the 
U2DIIF code: h/b=.05, and a,=.0873 radian, and the results were scaled 
appropriately. Scaling lift and moment results followed a linear relationship 
which implies that scaling the amplitude of motion by a factor results in scaling 
the forces up by the same factor. Drag follows a nonlinear relationship since 
drag is proportional to a*. Therefore, scaling the drag required multiplying the 
results from U2DIIF by the square of the scaling factor. After making all these 
conversions the U2DIIF harmonic output was compared to Bosch results in 
Figures 2.9 and 2.10. 

Excellent agreement is shown in plunge for both the drag coefficient, 
and the propulsive efficiency. This is not the case for pitch. Good agreement 
is shown for propulsive efficiency in pitch, but the drag coefficient does not 
agree well at the higher reduced frequencies. The reason for this is most likely 
the result of scaling errors. When small amplitudes are used for pitch motion 
the resulting drag coefficient magnitudes are very small ~.005. This is near 
the same magnitude as the accuracy of the code. Therefore, the inherent error 
in the code is scaled with the coefficient of drag magnitude. Plunge drag 
coefficient data are an order of magnitude larger than pitch and do not have the 


same scaling problem. 
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Figure 2.9 Average drag and efficiency for a NACAOOOS9 airfoil undergoing 
pure plunge oscillations at an amplitude of .05c. The results were scaled up 


to 1.0c (x100) for comparison to Bosch. 
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Figure 2.10 Drag and effictency for a NACAOOO9Y airfoil undergoing pure 
pitch oscillations at an amplitude of 5.7°. The results were scaled up to an 


amplitude of 57 ° (x100) for comparison to Bosch. 
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4. Power Extraction for Two Degrees of Freedom 

When an airfoil undergoes both pitch and plunge motion it is 
oscillating with two degrees of freedom. For this condition it is possible for a 
phase relationship to exist between the two types of motion which results in 
a positive energy input into the airfoil. When this phase relationship exists 
energy is extracted from the airstream and transferred to the airfoil. This is the 
basic mechanism behind aerodynamic flutter. This may also be used as an 
efficient source of mechanical energy as described by McKinney and Delaurier; 
The Wingmill: An Oscillating-Wing Windmill (Ref. 9]. The objective of this 
study was to compare the results obtained from the U2DIIF code for power 
extracted from airfoil motion in two degrees of freedom to the experimental 
results obtained in reference 9. First, the U2DIIF code was modified to perform 


the following harmonic motion, 
h=hsin(wf) , a = a sin(wf+d) (2.19) 
where, @ is the phase angle between plunging and pitching motions. Next, the 
equation for power extracted from the airstream is given by, 
P = h(Ncos« + (T,-D)sina) + aM (2.20) 
where, N = normal force, D = drag, T, = suction force, and M = moment 


about the pitch axis. Assuming small perturbations and noting the opposite 


signs of T, and D, we can write, 


25 


(7, - Dsina < Nceosa , Lift = N (2.21) 


Therefore, the equation for power is reduced to, 


P=Lh+ Mé (2.22) 


Equation 2.22 was used to obtain the power extracted during one 
cycle of motion. Experimental results are given in reference 9 for a wing 
undergoing the harmonic motion described in equation 2.19. The wing has a 
NACAOO12 airfoil section and a rectangular planform of 20cm by 105cm. In 
the experiment h, was kept constant at 6 cm, and q was set to 25°. The 
reduced frequency was given to be .361 based on the half chord and the phase 
angle was set to 90°. 

Figure 2.11 shows the airfoil motion for both angle-of-attack and 
translation as a function of time for one cycle. This plot shows the 90° phase 
difference between a and h. Figure 2.12 shows the aerodynamic forces in 
newtons and the power extracted in watts obtained from reference 9 and the 
U2DIIF code. Numerically, the experimental results differ from what we 
obtained by a factor of 2 both in the aerodynamic forces and the power 
extraction. They obtained a peak normal force coefficient equal to .897, while 
we obtained a peak lift coefficient of 2. This difference is due to the three- 
dimensional effects associated with the experiment combined with our 
approximation of the normal force coefficient. In reference 9 a lift-curve slope 


of C,, = 4.30/rad was used to correct the two-dimensional Theodorsen theory. 
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If we apply the same correction to the theoretical lift-curve slope of c,, = 27 
used in U2DIIF along with the normal force correction, we obtain a normal force 
coefficient equal to 1.23. This is slightly greater then the results of reference 
9. This is what we expected considering the low Reynolds number and added 
mass term used in the wind tunnel experiments which would tend to lower the 
results from theoretical predictions. 

The power extraction curve followed the same trend as the 
aerodynamic forces, and again was off by a factor of 2. After applying the 
same corrections as above we obtain a peak power equal to 16.6 watts 
compared to the experimental value of 13.5 watts. It was encouraging to see 
the predictions from U2DIIF produce the same shape and phasing as the 
experimental results. Furthermore, we accurately predict positive work (power) 


over most of the cycle, as is expected for a phase angle of 90°. 
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Figure 2.12 The aerodynamic forces (newtons) and power extracted (watts) 
for both the results obtained by McKinney [Ref. 9] and the results from the 
U2DIIF code respectively. 
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C. UNSTEADY PRESSURE DISTRIBUTIONS 


1. Introduction 

The purpose of this study was to establish relationships between 
steady state pressure distributions and unsteady pressure distributions for a 
NACAO0O12 airfoil using the panel code U2DIIF. The study centered on finding 
an unsteady pressure distribution on the upper surface of an airfoil undergoing 
a ramp motion which would match a pressure distribution from the same airfoil 
at some steady state value. Once the relationship was developed between 
steady and unsteady pressure distributions, the results could be used as an 
input in a direct boundary layer code. This would give a quick and efficient 
method of determining unsteady boundary layer profiles and the onset of 
dynamic stall using a simple panel code coupled with a direct boundary layer 


code. 


2. Theory 
The principal theory or hypothesis is that an airfoil undergoing a quasi- 
steady ramp motion can be analyzed as a steady-state process. The current 
thought process is that a ramp motion is a purely unsteady phenomenon which 
requires an unsteady analysis (ie. Navier Stokes) to determine the boundary 
layer and the onset of dynamic stall. In fact, as pointed out in reference 11, 
dynamic stall is strongly dependent upon airfoil geometry, Mach number, pitch 


rate, Reynolds number, state of the airfoil boundary layer, type of motion, etc. 
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Furthermore, the computation of unsteady, compressible flows can be a long 
and complicated process, even assuming the flow can be modeled 
appropriately. The next thought process Is to assume that the boundary layer 
behaves in a quasi-steady fashion. This implies that the boundary layer reacts 
nearly instantaneously to the pressure distribution. More precisely, the 
boundary layer will react to the pressure gradients. The stronger the adverse 
pressure gradient the more likely the boundary layer tendency is to separate. 
To check this theory, a study was done in reference 10 where the unsteady 
pressure distribution for a particular motion was input into a steady direct 
boundary layer code. The direct boundary layer code can then accurately 
predict the onset of separation or flow reversal, although the code breaks down 
after separation due to the formulation of the direct boundary layer problem 
[Ref. 10]. The onset of separation for different pitch rates was compared to 
experimental results obtained by Chandrasekhara, Carr, and Ahmed at NASA 
Ames Research Center [Ref. 11]. This process covers all the variables stated 
above that influence dynamic stall with some limitations. First, the airfoil 
geometry, type of motion, and pitch rate are covered in the unsteady panel 
code. The effects of Reynolds number, and Mach number are covered in the 
direct boundary layer code because this code was modified to include the 
Prandtl-Glauert compressibility correction factor. The primary limitations in this 


analyis are that the three-dimensional effects are not considered, and at higher 
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pitch rates the flow becomes transonic making the compressibility correction 
less accurate. 

Once it has been established that the boundary layer behaves in a 
quasi-steady fashion for certain flows, a relationship can be obtained that links 
the quasi-steady pressure distribution on the upper surface to a purely steady 
state pressure distribution. For this to be true we must assume that the entire 
flow field is quasi-steady. This requires keeping the ramp motion sufficiently 
slow to satisfy this assumption. When an airfoil moves in a quasi-steady ramp 
motion the pressure field will require a finite amount of time to react as 
explained in the first section. For this reason the pressure distribution on the 
upper surface of the ramping airfoil will match a steady state pressure 
distribution at some earlier AOA. The purpose of this study was to verify this 
hypothesis and to determine the relationship between steady and quasi-steady 


flow fields. 


3. Description 
This study involved using the U2DIIF code for NACAOQO12 airfoils 
undergoing a ramp motion at various rise times to determine the relationship 
between steady and quasi-steady flows. It was first necessary to determine 
the range of rise times for which the assumption of quasi-steady flow would 
hold true. The rise time is defined as the time required for the airfoil to ramp 


from the initial AOA to the final AOA. The rise time is expressed in 
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nondimensional units by Teng [Ref 2.] and is defined as the time required for 


a fluid particle to travel from the leading edge to the trailing edge, 


f= (2.23) 





Since experimental results can normally be assumed quasi-steady due to 
mechanical limitations, the experimental results in reference 11 were used as 
a guide. In these experiments the ramp motion goes from O to 5” in 
0.00167sec. The fluid particle takes 0.0005 sec to go from the leading edge 
to the trailing edge. Therefore, the nondimensional rise time is equal to 3.34. 
This is equivalent to a nondimensional pitch rate as defined in reference 11 of 
.027 units. Since this study involved ramp motions from O to 15”, rise times 
on the order of 10 units were considered quasi-steady. 

The study involved running the U2DIIF code for a particular rise time. 
The pressure distribution output was then compared to the steady state 
pressure distribution for AOA‘s of 6, 8, and 10 degrees. The comparison 
continued until the best match was found between the unsteady pressure 
distribution and the steady state distribution. In order to compare the pressure 
gradients on the upper surface, the pressure distribution from the U2DIIF output 


was differentiated using a four point central difference scheme. 


4. Results 
Figures 2.13 and 2.14 show a typical comparison study of the 


pressure distributions and pressure gradients. The steady state distribution is 
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plotted for a particular AOA along with several unsteady distributions for 
different AOA’‘s until a ‘best’ match is found. In this particular study the match 
fora 10° AOA steady distribution was found at the unsteady AOA of 13.34° 
for a relatively fast rise time of 3 units. Figure 2.15 shows the AOA and the 
aerodynamic forces as a function of time for this ramp motion. This important 
result verifies that there does exist an unsteady pressure distribution that will 
match a steady state pressure distribution on the upper surface. 

One of the first results discovered in this study was that when the 
rise times were decreased from 10 to 3 units the best match became 
increasingly off. This was to be expected as the ramp motion increased to 
where the quasi-steady assumption would be invalid. It was found that rise 
times faster than 3 units would not match steady state pressure distributions. 
To show the trend, Figures 2.16 and 2.17 show the best match of a 10° AOA 
steady pressure distribution and the corresponding unsteady pressure 
distribution for a rise time of 10 and 3 units respectively. As can be seen, the 
match becomes progressively worse with the faster rise times. 

The numerical results for this study are tabulated in Table 2.1. 
Certain trends can be observed from these results. First, the steady state lift 
coefficient at a particular AOA is always greater than the corresponding 
unsteady lift coefficient at the same AOA. This is consistent with Theodorsen 
theory. Next, the peak unsteady lift coefficient is always greater than its 


corresponding steady state value and increases with increasing pitch rate. This 
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is consistent with the increased maximum AOA and corresponding lift obtained 
from dynamic motion over steady state results as shown in reference 11. 

The AOA for which a particular match was found also followed 
certain trends. First, the delta AOA at the smaller steady state angle-of-attack 
(6°) was almost independent of the rise time. This is equivalent to stating that 
it is independent of the pitch rate. A nominal value of two degrees was found 
between the unsteady and steady distributions for the smaller steady state 
AOA‘s. Next, the higher the steady state AOA, the more dependency on pitch 
rate was found. At 10° AOA steady the delta AOA went from 2.13 degrees 
for a rise time of 10 units to a delta AOA of 3.2 degrees at a rise time of 3 
units. 

When comparing this data to the experimental results of reference 11, 
it was not straight forward. The best we could do was compare the effects 
of pitch rate on the onset of separation. Reference 11 notes that when the 
pitch rate is doubled from .03 to .05 units the onset in separation is delayed by 
2° (from 14° to 16°). Although we do not compute separation, we can 
compare the delay in matching the steady state values at an AOA of 14° for 
the pitch rates used in reference 11. The results of this study are tabulated at 
the end of Table 2.1. Ascan be seen, we achieve the same delay of 2° when 
we double the pitch rate for a rise time of 14 to 7 units with a ramp change of 


20°. 
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TABLE 2.1 Results of unsteady to steady pressure distribution matching. 


Steady Steady Unsteady 
AOA Cc. AOA AOA 

Rise Time = 10 units (O° - 15°) 

835 

960 | tora | 214 | 1.067 

10 1.254 
Rise Time = 8 units (O° - 15°) 

723 

1.07 

1 1.26 


Rise Time = 6 units (O° - 15°) 


902 


| 


1.104 
10 1.281 
Rise Time = 4 units (O° - 15°) 
| 960 | 10.36 | 236 | 1.18 
= — a 

Rise Time = 3 units (O° - 15°) 
1.18 
960 | toase | 2.43 | 1.26 
10 1.35 


Rise Time = 14 units (O° - 20°) 


14 


Rise Time = 7 units (0° - 20°) 
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Figure 2.15 Airfoil motion and aerodynamic forces for a NACAOO12 airfoil 
undergoing a ramp motion from O to 15° in arise time of 3 units. 
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lil. TWO AIRFOIL ANALYSIS 


A. COMPUTER CODE USPOTF2 
The computer code written by Pang [Ref. 12] was the primary method 
used in the following studies of two airfoils undergoing harmonic motion. The 
computer code, USPOTF2, was written to solve the potential flow for two 
airfoils executing unsteady motions in an inviscid, incompressible flow medium. 
This code is an extension of the single airfoil code, U2DIIF [Ref. 2], which uses 
the technique known as Pane/ Methods for steady flow and extends it to 
unsteady flow by introducing a wake model. This creates a non-linear problem 
due to the continuous shedding of vortices into the trailing wake. Furthermore, 
the presence of the second airfoil introduces a set of non-linear coupled 
equations for the Kutta condition. The solution requires an iterative procedure 
to compute the two vorticity strengths. Although no attempt is made here to 
reproduce the work by Pang, a general list of the modifications required to 
enhance the original code to a two-airfoil code is shown here. This will provide 
the reader with the necessary information to understand the modifications 
made by this author: 
1. The establishment of five frames of reference: one fixed inertia frame of 


reference (global), two moving local frames of reference and two frozen 
local frames of reference. 
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2. Reformulation of the two Kutta conditions which are coupled non-linear. 


3. The creation of a new subroutine (NEWPOS) which transforms all 
coordinates in either of the two respective local frames of reference to 
the global frame of reference. 


4. The introduction of a more accurate method to obtain the velocity 
potential by integrating the velocity over smaller panels on the airfoil. 


5. Extension of the influence coefficient to include the effects of the 
second airfoil with its own peculiar wake. This also requires an 
introduction of an additional influence coefficient, that on the wake 
element due to the wake element from the other airfoil. 


6. The program is restricted to the following types of motion: 


In-phase and out-of-phase Step Input 


In-phase and out-of-phase Modified Ramp Input 


In-phase and out-of-phase Translational Harmonic Oscillation 


@ in-phase and out-of-phase Rotational Harmonic Oscillation 


Sharp Edge Gust Field Penetration 
1. Modifications (USPOTF2A) 


a. Program Output 
The code was originally written on the IBM mainframe at the 
Naval Postgraduate School. It has been transferred to the Stardent machine of 
the Department of Aeronautics and Astronautics. This allowed for increased 
storage Capacity and a decrease in run times. For this reason, the total number 
of panels used to describe both airfoils has been crasiiad from 200 to 400. 
The program was also transferred to the Iris workstation in the department. 


Although the program would compile on the workstation, this author could not 
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get the program to run on the Iris. It was believed to be a memory storage 
problem, but a thorough investigation was not completed. This program 
produces an extremely large amount of output to the screen. It is usually 
convenient on a Unix based machine to redirect screen output to a file during 
program runs. This file can get very large for long program runs, and may 
cause a problem in the temporary storage capacity. For this reason, the logical 
variable ‘output’ was added to the input file, FOROO1.DAT. When this variable 
is set to fa/se, most of the screen output will not be printed, thus reducing the 
size of the output file. Furthermore, the important output has been separated 
into different files for easier analysis, and is not affected by the variable output. 


The following list describes the input/output files and the data they contain: 


ie FOROO1.DAT: This is the input file (see Figure 3.1). 

De FOROO2.DAT: This is for user supplied airfoil coordinates, if 
desired. 

3. FOROO3.DAT: This file contains the global coordinates of the first 


airfoil at each time step. 


4. FOROO4.DAT: This file contains the global coordinates of the 
second airfoil at each time step. 


D: FOROO?7.DAT: This file contains the lift, drag, and moment 
coefficients for both airfoils at each time step. 


6. FOROO8.DAT: This file contains the pressure coefficients for the 
first airfoil at each time step. 


ye; FOROOY.DAT: This file contains the pressure coefficients for the 
second airfoil at each time step. 


at 


10. 


11. 


EZ. 


13. 


14. 


FORO10.DAT: This file contains the first airfoil’s core vortex (wake) 
positions at each time step. 


FORO11.DAT: This file contains the second airfoil’s core vortex 
(wake) positions at each time step. 


FORO12.DAT: This file contains the AOA at each time step. 


FORO13.DAT: This file contains the DHY (translational motion) at 
each time step. 


FORO14.DAT: This file contains the required input for the Phase 
program. 


FORO20.DAT: This file contains the computed average pressure 
coefficient for the first airfoil. 


FORO21.DAT: This file contains the computed average pressure 
coefficient for the second airfoil. 


Program Corrections 


There were several errors noted in the original program and 


corrected. First, the code cannot be executed with only one airfoil as implied 


by the input variable NAIRF, or number of airfoils. Two airfoils must be defined 


to run the program. If single airfoil results are desired, one must position the 


second airfoil at a sufficient distance (approx. 30 chord lengths) away from the 


first airfoil to ensure no interference effects. This procedure was used to 


compare the USPOTF2 code against the single airfoil code, U2DIIF. At first, 


the results did not match. This was when an error in the subroutine Press was 


discovered and corrected. Now, if one selects NGIES equal to 0.0, which is the 
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unsteady Kutta condition of equal pressures at the trailing edge panels, the two 
codes produce the same results. 

Next, an error in the convergence criteria was detected. 
Originally, the code would not reset the variable TOL, which is used by the 
code to set up the tolerance criteria for convergence. If, for a particular time 
step, the code required a higher tolerance to converge, the code would use this 
higher tolerance for follow-on time steps. This was corrected by resetting the 
variable TOL to its original value at the end of each time step. 

The original code computes the moment coefficient about the 
leading edge. Although this is not an error, it is usually desirable to have the 
moment coefficient about the pitch axis. This change was made in the 


modified code. 


c. Airfoil Motion 

This was the largest modification made to this code. As stated 
earlier, the original code was restricted to in-phase and out-of-phase motion. 
The code was modified so the airfoils can move independently of each other. 
This was done by adding the following variables into the code: DALP2, TCON2, 
FREQ2, and PIVOT2. The entire logic for airfoil motion was rewritten to include 
the new variables which allowed the airfoils to move independently. 
Corrections for this modification were made as needed in the subroutines. 
Figure 3.1 shows a sample input file, FOROO1.DAT, required for the modified 


version of USPOTF2. Figures 3.2 and 3.3 demonstrate some of the new 
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capabilities of the modified code. In Figure 3.2, the first airfoil is stationary at 
10° AOA, while the second airfoil is oscillating in pure translation (DHY =.1) 
at a reduced frequency of 2.0. The aerodynamic forces clearly show an 
influence of the second airfoil’s oscillation (drag coefficient) on the first airfoil’s 
lift coeffiecient. In Figure 3.3, both airfoils are oscillating in pure pitch (DALP 
= 10°), but at different reduced frequencies. The top airfoil is oscillating at 
a reduced frequency of 4.0, while the bottom airfoil is oscillating at a reduced 
frequency of 1.0. The aerodynamic forces are shown for 4 cycles of the top 
airfoil which is equivelant to 1 cycle of the bottom airfoil. The interaction of 
both airfoils can be seen in the aerodynamic forces with the stronger influence 
from the fast oscillating top airfoil. 

Another modification was made in the original code to allow for 
different size airfoils to be analyzed at the same time. The original code 
nondimensionalizes both airfoil’s chord length to 1.0. The code was modified 
to include a new variable SCALE. This variable scales the second airfoil by the 
amount specified as a percent of the first airfoil. Figure 3.4 demonstrates the 


capability of this added feature. 
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Aug 17 1993 16:37:30 stdin 


1 NUMBER OF LINES FOR TITLE 
2 1 
3 TWO NACAON012 OCSILLATING AIRFOILS 
4 
5 IFLAG NLOWER NUPPER 
6 00, 30, 30 
Z NATREOG, XSHIET eYSHIET, SCALE 
8 Z; 12° sO 1 8, 
9 NACA AIRFOIL TYPE, 
10 LZ 
il V2 
#2 ALF], ALE2Z, DALP1, DAEP2, TCONT.. FeONnZ. 
2 O20, U0, “Poe, 6.0 eel or ape ap 
14 FREQ], FREQ2, PIVOT1, PIVOT2 
1s 020-5 .ac0 0s O..0% (OFF Os 
16 UGUST, VGUST, DELHX1], DELHX?, DE! HY], DELHY2, PHASE], PHASE? 
17 0.05" 0.0. 0 2G. O70: O20; On2, O20S 207 
18 1G ae DTS, DTS2, TOL, TADJ, SCL, SCM, SGAM, NGIFS 
19 bed, 005; 7020. 22009) 0 0 O0F 0 0e O20; 0 
20 STEADY OUTPUT 
| 21 false, true 
22 
23 STEADY--TRUE IF ONLY STEADY SOLUTION. FALSE OTHERWISE. 
24 OUTPUT--TRUE IF YOU WANT COMPLETE OUTPUT TO SCREEN. 
25 IFLAG : O IF AIRFOIL IS NACA XXXX OR 730xx 
26 1 OTHERWISE. 
2? NLOWER : NO. OF PANELS USED ON BOTH AIRFOIL LOWER SURFACES. 
28 NUPPER : NO. OF PANELS USED ON BOTH AIRFCIL UPPFR S!'RFATES, 
29 NAIRFO : NUMBER OF AIRFOILS. 
30 XSHIFT : RELATIVE X DIST. FROM 2 AIRFOIL PIVOT FOSITION &RT 
31 GLOBAL COORDINATE SYSTEM. 
SZ YSHIFT =: RELATIVE Y DIST. FROM 2 AIRFOIL PIVOT POSITICN eer 
33 GLOBAL COORDINATED SYSTEM. 
34 NACA AIRFOIL TYPF * ENTER NACK 4 OR-5 DIGIT COpE FOR AlRec tus 
35 IF NCT A NACA AIRFOIL, SUPPLY AIRFOYTI. 
36 X{I),Y(1) COORDS. FOR BOTH AIRFCILS IN 
37 PITESCODE s2, 
38 ALP1/2 : INITIA!. ANGLE OF ATTACK FOR AIRFOILS IN DEGREES. 
} 39 DALP1/2 CHANGE IN AOA IN DEGREE FOR NON OSCILL. MOTION. 
40 MAX AMPLITUDE OF AOA IN DEGREE FOR ROT. HARMONI’” MCTION. 
41 TCON]/2 : NON-DIMENSIONAL RISE TIME (Vinf.t/C) OF ACA £OR 
42 MODIFIED RAMF CHANGE IN AOA, 
43 FREQI1/2 : NON DIMENSIONAL OSCILL. (WC/Vinf.) FOR HARMONIC MOTIONS. 
44 PIVOT1/2: LENGTH FROM LEADING EDGE TO PIVOT POINT FOR LOCAL SYSTEM. 
4s (THE GLOBAL SYSTFM’S ORGIN IS THE FIRST AIRFOILS PIVOT POSITION) 
46 UGUST : MAG. OF NON-DIM. GUST VELOCITY ALONG GLOBAL X DIRECTION. 
47 VGUST : MAG. OF NON-DIM. GUST VELOCITY ALONG GLOBAL Y NOIRFCTION. 
48 DELHX1/2: NON-DIM. TRANSLATICNAL CHCORDWISE AMPLITUDE. 
| 49 DELHY1/2: NON-DIM. TRANSLATIONAL TRANSVERSE AMPLITUDE (plunging). 
S¢ PHASE1/2: PHASE ANGLE IN DEGREE BETWFEN CHORDWISE AND TRANSVERSE 
31 TRANSLATIONAL OSCILIL. WITH THE LATER PER Os THAT Al RibnOiles 
52 TE : FINAL NON-DIM. TIME TO TERMINATE UNSTEADY FLOW SOLUTION. 
53 DES ¢ STARTING TIME STEP FOR NON-OCIILL MOTIONS(TADJ=C). 
54 NO. OF COMF'JTAICNAL STEPS PER CYCLE FOR HARMONIC MOTION 
ate: (FOR 2 FREO OCILL. IT USES THE LARGEST FEREC) 
56 BASELINE TIME STEP FOR ALL MOTIONS(TADJ NOT =0) 
5? DTS2 : STARTING NON-DIM TIME FOR SECOND AIRFOIL MOTION TO BFGIN. 
58 (O TO BEGIN MOTION AT THE SAME TIME). 
$9 TOL : TOLERANCE CRIERION FOR CONVERGENCE FOR (Uw) and (Vw)k. 
60 PADS ¢ FACTOR BY WHICH DIS WILL BE ADJUSTED. 
61 SCLA ¢ STEADY LIFT COEFF. FOR THE SINGLE ATRFCIL Al Tit Sees or 
62 SCM : STEADY MOMENT CCEFF. FOR THE SINGLE AIRFOIL. 
63 SGAM : STEADY VORTICITY STRENGTH FOR THE SINGLE AIRFOIL. 
64 NGIES OPTION TO CHANGF THF UNSTEADY KUTTA CCNDITION. 
65 O EQUAL PRESSURE AT THE TRAILING EDGE PANELS. 
66 1 EQUAL TANGENTAL VELOCITIES AT THE TRAILING ENGE PANELS. 


Figure 3.1 Sample input file, FOROO1.DAT, required for the modified version 


on USPONEZ: 
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Page 1 


Wake Pattern 


“bewveacaee™ 





mie 


Aerodynamic Forces 


Magnitude 





0 2 4 6 
NONDIM. TIME 


Figure 3.2 Wake pattern and aerodynamic forces. AF1-10 deg steady, AF2-0 
deg unsteady (DHY =.1,k = 2) 
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Figure 3.3 Wake pattern and aerodynamic forces. AFl-unsteady (DALP =10 
deg, k=1), AF2-unsteady (DALP = 10 deg, k=4). 
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d. Limitations 

There are several important limitations that must be understood 
thoroughly before attempting to use this program. First, the core vortices that 
make up the airfoil’s wake must not intersect the panels that define the airfoil 
geometry. This situation is possible if one places the second airfoil in the wake 
of the first airfoil. The singularity nature of a core vortex intersecting a source 
can sometimes cause erroneous output. This can be noted in two ways. One 
way is to watch the screen output during the program run. If this situation 
arises, the code will not converge until the tolerance has been changed several 
times and is easily observed. Another way to identify this problem is to 
examine both the force, and moment data along with the wake patterns at the 
conclusion of the run. If the data or wake patterns are discontinuous, then 
there is a problem. This problem can be less severe if the panel density is 
increased and the time step is decreased. Unfortunately, if the wake Is strong 
the interference effects are too great and the solution is suspect. Figure 3.5 
shows the wake pattern and aerodynamic forces of a stationary airfoil in the 
strong wake of an oscillating airfoil. This example clearly points out the biggest 
limitation with this program. 

Other limitations include a sensitivity to panel density, time step, 
and amplitude of motion. When adjusting these parameters one must watch 


the screen output to ensure the program can converge at each time step. 
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Figure 3.5 Wake pattern and Aerodynamic forces showing vortex interference 
effects. AF1 unsteady, AF2 steady. 
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B. COMPARISON TO EXISTING COMPUTATIONAL RESULTS 

D.D. Liu and Z.X. Yao published results for vortex/wake flow studies in 
reference 13. For these results a numerical scheme was developed at Arizona 
State University that uses Panel Methods combined with an unsteady model 
that sheds vortices into the wake at each time step. The primary difference 
between this computer code and USPOTF2A is the treatment of the unsteady 
Kutta condition. Liu and Yao use a linearized model following the approach by 
Kim and Mook, where USPOTF2A uses the non-linear model developed by Basu 
and Hancock. 

Two case studies were chosen to compare the results. First, consider the 
case where two NACAOQO12 airfoils undergo a step change in AOA from O to 
5°. The second airfoil is placed .5c behind and .2c above the first. The results 
for wake patterns and aerodynamic forces for USPOTF2A are shown in Figure 
3.6. Liu and Yao have revised their results from reference 13. The current 
results are published here in Figure 3.7. The next case study involved the first 
airfoil undergoing a step change in AOA from O to 5” while the second airfoil 
remained stationary. Again, the second airfoil was placed .5c behind and .2c 
above the first. The results obtained by USPOTF2A are shown in Figure 3.8. 
The revised results obtained by Liu and Yao are shown in Figure 3.9. 

It is clear from the results that both codes follow the same trends. They 


do not produce the exact same results, as expected, due to the different 
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formulation of the problem. Both results essentially produce the same 
interference effect and approach steady state conditions. This example is a 


good verification for both numerical codes. 
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Figure 3.6 USPOTF2A wake pattern and lift history output for two 
NACA00O12 undergoing a step change in AOA from O to 5°. 
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Figure 3.7 Liu’s output of lift histories for two NACAO012 airfoils undergoing 
a step input from O to 5°. 
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Figure 3.8 USPOTF2A wake pattern and lift history output for two 
NACAO012 airfoils. AF-1 undergoing a step change in AOA from O to 5°. 
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Figure 3.9 Liu’s output of lift histories for two NACAO0O12 airfoils, with the 
fist undergoing a step input from O to 5°. 
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C. PROPULSIVE EFFICIENCY 


1. Introduction 

The purpose of this section is to develop relations for propulsive 
forces and efficiencies for two airfoils in tandem. There are two configurations 
discussed in this section (Fig. 3.10). One, the first airfoil oscillates while the 
second remains stationary. Two, the first airfoil remains stationary, while the 
second airfoil oscillates. The computational study will only examine the second 
configuration where the airfoil oscillates in either pitch or plunge. Due to the 
limitations of the code USPOTF2A, the case where the first airfoil oscillates is 
discussed with limited computations. The results of the numerical output were 


compared to the analytical results obtained by Bosch in reference 8. 


Configuration _1 Configuration 2 


— a — aa 


Plunge Motion Plunge Motion 


— — oa 


Pitch Motion Pitch Motion 





Figure 3.10 The two configurations studied. 
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2. Theory 

The theory for the development of propulsive forces and efficiencies 
for two airfoils in tandem follows closely to the single airfoil case presented in 
Chapter Il. The relationships for work done to the fluid by the two airfoils 
remain the same as for the single airfoil case since only one of the airfoils is in 
motion. The relationships for work done by plunge motion is given by equation 
2.15 and for pitch motion by equation 2.17. The primary difference between 
the single airfoil analysis and the two airfoil analysis is in the relationship for 
efficiency. Although only one airfoil is in motion, both airfoils are capable of 
producing thrust, and the expression for efficiency becomes: 


_ (D,+D,)-v 
W 


n (3.1) 


where D,, and D, are the drag or propulsive forces from the first and second 
airfoil, respectively. 

From the expression for efficiency in equation 3.1, the differences 
between the two configurations can be examined. In the first configuration, 
the pressure disturbances from the first airfoil are felt by the second airfoil 
downstream, or in the first airfoil’s wake. This will produce a much stronger 
effect than the second configuration; where the pressure disturbances must 
travel upstream to have an effect on the stationary airfoil. Since the production 


of thrust by the stationary airfoil is entirely dependent on the influence of the 
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oscillating airfoil, it is clear that the first configuration will have a greater 
efficiency than the second configuration. Furthermore, as stated in chapter 2, 
plunge motion has a stronger influence than pitch motion; and the efficiencies 
associated with plunge motion are greater than those associated with pitch. 
Bosch concluded this in reference 8. Efficiencies for the first configuration near 
.9 were obtained for both types of motion with reduced frequencies between 
1 and 2. The second configuration reached an efficiency near .5 for plunge 


motion and .1 for pitch motion for the same range of reduced frequencies. 


3. Comparison to Flat Plate Theory 

Analytical results for two flat plates undergoing pure pitch and pure 
plunge motion for both configurations were computed by Bosch [Ref 8.]. To 
compare the results using USPOTF2A, the same conversions and scaling of the 
motion as for the single airfoil case were applied. Special consideration was 
given to the panel density and the time step to reduce the effects of the first 
airfoil'Ss wake impinging upon the second airfoil. It was shown that when the 
first airfoil’s wake is ‘weak’, it has a minimal effect on the solution or wake 
pattern of the second airfoil. The aerodynamic output was carefully observed 
proceeding each case study to ensure a smooth and continuous curve; this 
would confirm the impinging wake did not adversely effect the results. The 
time dependent output was then converted to harmonic output following similar 
procedures for the single airfoil case using the two airfoil phase-shift code 


presented in Appendex B. Then, the magnitudes of lift and moment were 
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divided into real and imaginary parts and plotted against reduced frequency. 
The results for pitch motion are shown in Figure 3.11; the first plot (C, ’) 
represents the real part, and the second plot (C, ’’) represents the imaginary 
part. The agreement is particularly good at the low frequencies and begins to 
separate at the higher frequencies. 

The comparison for thrust and efficiency was done using the second 
configuration with plunge motion. The expressions for efficiency are given by 


the following relations: 


— Cp + Cp, 2(Cp + Cp) 


= yt De z (3.2) 
eee Kim(Cuy "punge ~~ kim(C,) 


Figure 3.12 shows the propulsive force and efficiency comparisons. Here, the 
results match well over the entire frequency range with the computed 
efficiency always being less than that predicted by Bosch. This was 
anticipated because the non-linear effects in the USPOTF2A code would result 
in a lower efficiency. Furthermore, the consequence of plunge motion 
producing higher aerodynamic forces for the same frequency as pitch motion, 
results in less scaling errors. 

A final study was done to try and achieve the results by Bosch for the 
first configuration where the efficiencies are much higher. Unfortunately, as 
Stated this code is unable to converge when the discrete vortices collide with 


the airfoil. Therefore, the second airfoil was placed O.2c in front and 1c above 
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the stationary airfoil. In this configuration, the trailing vortices from the 
oscillating airfoil would not interfere with the stationary airfoil. The results 
were recorded for plunging motion (DHY = .05) at areduced frequency of 4.0, 
and then the oscillating airfoil was brought incrementally closer to the 
stationary airfoil. This continued until the vortices interacted at 0.3c above the 
Stationary airfoil. The efficiency went from 7=.48, to 7=.6 at 1c and .3c 
above the stationary airfoil, respectively. This study shows the correct trend 
for the efficiency under the limitations of the code. 

Overall, the USPOTF2A program was shown to be an effective way 
to predict the aerodynamic forces, and efficiencies for two interacting airfoils. 
Unfortunately, the limitations of the program prevent some interesting case 
studies of strong wake interaction. Finally, one must be cautioned to observe 


the output with care when studying airfoils in close proximity to the wakes. 
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Figure 3.11 Lift coefficient as a function of the reduced frequency k. 
Pitching motion using configuration 2. 
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Figure 3.12 Drag coefficient and propulsive efficiency as a function of the 
reduced frequency k. Plunging motion using configuration 2. 
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D. UPWIND INFLUENCE STUDY 


1. Introduction 

It has been of interest to understand and study the effects of an 
oscillating airfoil in the vicinity of a stationary airfoil. The purpose of such a 
study is to see if there exists a configuration that would enhance the lift on the 
Stationary airfoil. This is of primary interest in the post-stall region, where an 
oscillating airfoil may have the beneficial effect of delaying boundary layer 
separation and increasing the steady state lift. In theory, the oscillating airfoil 
will produce a certain amount of thrust, thereby increasing the flow velocity 
and creating a more favorable pressure gradient on the stationary airfoil. 
Although this study is done with the inviscid code, USPOTF2A, useful 
information on the pressure distributions and gradients are obtained from such 
an analysis. This information can be used to understand the degree of 
influence and trends in the integrated forces that determine the effects on a 
certain configuration. Due to the limitations of the code, only configurations 


where the oscillating airfoil was behind the stationary airfoil were considered. 


2. Description 
This study used the modified code, USPOTF2A, to determine the 
effects of an oscillating airfoil on a stationary airfoil at a high AOA. Two 
NACAOOOS airfoils were used for this study with the first remaining stationary 


at 10° AOA. The second airfoil was set to oscillate in plunge (DHY =.2) ata 
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reduced frequency of 4.0. The amplitude of motion, DHY, was limited to .2 
because of the wake interference problem. The high reduced frequency was 
chosen to get the maximum effect with the relatively small amplitude of motion 
used in this study. The study involved changing both the position and mean 
AOA of the oscillating airfoil; then recording the change in pressure distribution 
and average lift of the stationary airfoil over steady state values. Two positions 
for the oscillating airfoil were studied. First, the oscillating airfoil was placed 
slightly behind (XSHIFT = .1c) and above (YSHIFT = .6c) the stationary airfoil. 
Here, the mean AOA of the oscillating airfoil was changed from O° to 10°. 
The other position considered was placing the oscillating airfoil slightly behind 
(XSHIFT = .1c) and below (YSHIFT = -.6c) the stationary airfoil. Here, the 


mean AOA of the oscillating airfoil varied from O°, 5°, 10°, and 15°. 


3. Results 

The numerical results of this study are tabulated in table 3.1. In the 
first configuration, where the oscillating airfoil was located above the stationary 
airfoil, the influence was minimal. The increase in lift ranged from 7% to 10%. 
The important result was the influence of AOA. For a greater mean AOA on 
the oscillating airfoil, the delta lift on the oscillating airfoil was increased 
dramatically along with the influence on the stationary airfoil. Furthermore, the 
production of thrust decreased with increasing mean AOA. The results of the 
second configuration, where the oscillating airfoil was below the stationary 


airfoil, proved to be more effective. The increase in lift ranged from 10% to 


68 


20%. Again, the influence of mean AOA followed the same trends as the first 
configuration. 

A typical case study is displayed in Figures 3.13 and 3.14. Figure 
3.13 shows the steady state pressure distribution super-imposed on the 
average pressure distribution for one cycle of motion. Furthermore, the 
pressure gradients for both cases are displayed on the same graph. Figure 3.14 
shows the aerodynamic forces for two cycles of motion. The second cycle 
was analyzed for average data. It is interesting to note the cyclic behavior of 
the stationary airfoil’s lift coefficient which is designated as AF-1. Clearly, the 


influence of the oscillating airfoil on the stationary airfoil is shown in this graph. 


4. Conclusion 

Several important conclusions can be made from the above results. 
First, there is certainly a favorable influence from an oscillating airfoil in the 
vicinity of a stationary airfoil. Furthermore, the AOA and position of the 
oscillating airfoil have a strong influence on the effectiveness of the 
configuration. It was shown that the most effective configuration was having 
the oscillator at a high mean AOA and below the stationary airfoil. Finally, it 
can be concluded that the upwind influence can only produce a maximum of 
17% change in lift for an unstalled mean AOA of 10°. Unfortunately, the 


limitations of the program prevented studying many interesting configurations. 
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TABLE 3.1: Numerical results of oscillating airfoil study. 
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IV. PROPULSIVE FLOW VISUALIZATION EXPERIMENT 


A. INTRODUCTION 

The purpose of this experiment was to document the production of thrust 
by the wave propeller originally built at the Naval Postgraduate School by Carl 
Dane [Ref. 15]. This was a preliminary experiment to better understand the 
vortex pattern produced by a wave propeller, and to examine the production of 
thrust using smoke flow visualization techniques. 

An explanation of what constitutes a propulsive vortical signature along 
with smoke flow visualization of the propulsive vortical patterns is given in 
reference 16. In this reference, the explanation is given by contrasting the 
vortical pattern produced by a cylinder (drag) with the vortical pattern produced 
by a pitching airfoil (thrust). The cylinder produced a vortical street where the 
top row of vortices rotated clockwise and the bottom row of vortices rotated 
counterclockwise. This pattern induces a velocity component in the upstream 
direction (Biot-Savart law). In contrast, the pitching airfoil produced a 
counterclockwise rotating vortex street on the top row and a clockwise rotating 
vortex street on the bottom row. This pattern induces a velocity component 
in the downstream direction. Reproduction of the flow visualization data from 


reference 16 is shown in Figure 4.1. 
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Figure 4.1 Results from reference 16. Top vortex street shed from a circular 
cylinder (drag). Bottom vortex street generated by an airfoil pitching about the 
quarter-chord (thrust). 


B. THEORY 

A comparison was done using the incompressible panel code, U2DIIF. The 
purpose of this study was to examine the vortical pattern produced by the 
panel code, and determine if the vortical signature matched experimental 
results. The input to the panel code was set up to best match the conditions 
of the experiment described in the next section. The panel code was run using 
a plunge amplitude, h,/c, equal to .364, a reduced frequency of 2.5 and a0O> 
mean AOA. The results of the vortical pattern are shown in Figure 4.2. Aside 
from the starting vortex, this is clearly a thrust producing vortical street. 


Furthermore, the vortical pattern is similar to that produced by the experiment 
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shown in Figures 4.10 and 4.11. Finally, the theory predicts a drag coefficient 


(thrust) of -.376. 


C. EXPERIMENTAL SETUP 


1. Wave Propeller 

The wave propeller, originally proposed by Wilhelm Schmidt [Ref. 14], 
used in this paper was a modification of the original construction made by Car! 
Dane [Ref. 15]. The original mechanism was designed to perform a circular 
motion while holding the wing at a constant angle-of-attack. The primary 
modification was the construction of a new wing. The new wing was made 
from a NACAOQOQ12 airfoil section and consisted of a 5.5’’ chord and a 22” 
span. The wing was built from a foam core covered with a thin plywood skin 
and finished with a layer of glass-epoxy composite for added fatigue strength. 
The wave propeller’s drive mechanism was rebuilt using added bearings and 
tighter fittings to allow for less binding at high rotational speeds. The final 
configuration for the wave propeller is shown in Figure 4.3. The drive motor 
used was a reversible DC motor rated at 24 volts, and 5 amps. The power 
supply used was regulated voltage DC power supply rated at 0-35 volts, and 
2.5 amps. The motor and power supply combination could provide the wave 
propeller with a maximum rotational speed of 1500 revolutions per minute 


(rom) when placed horizontal. 
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2. Wind Tunnel 

The wind tunnel used in this experiment was a very low speed, low 
turbulence smoke tunnel at the Naval Postgraduate School. This smoke tunnel 
is of indraft type and was designed to be a scaled model of an existing smoke 
tunnel at the NASA Ames Research Center. It is made of plexiglass walls and 
has a contraction ratio of 2.8:1. The motor provides wind tunnel velocities 
between O and 10 feet per second (fps). The smoke was created using a 
Rosco smoke generator and piped into the tunnel as a single stream for the 
streak line flow visualization experiment. The smoke was directed into a rake 
for the flow field visualization experiment. Figure 4.4 is a photograph of the 


wind tunnel and smoke rake used in this experiment. 
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Figure 4.3. Wave Propeller. 





Low speed smoke tunnel. 


Figure 4.4 
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D. TEST PROCEDURE 

Testing was done in the low speed smoke tunnel for two different 
conditions. First, with the wind tunnel off (O velocity), and second with the 
wind tunnel operating at 10 fps. In the first configuration with the wind tunnel 
off, the wave propeller was placed in the wind tunnel and the tunnel was filled 
with smoke; then, the wave propeller was turned on to its maximum rotational 
speed of 1500 rpm and a rotation diameter of 2 inches. The purpose was to 
see if the wave propeller would draw the smoke through the tunnel like a fan, 
thus showing the production of thrust by the wave propeller. 

In the second tunnel condition, the smoke tunnel was turned on to its 
maximum power bringing the tunnel speed to 10 fps (chord Reynolds number 
of 26,500). The wave propeller was oscillated at three different rotational 
speeds: 165 rpm, 620 rpm, and 1085 rpm. This equated to reduced 
frequencies of ./92, 2.97, and 5.2 respectively. Furthermore, the angle-of- 
attack of the wave propeller was varied between O, 5, 10, and 20 degrees with 
a rotation diameter of 2 inches. For these configurations, smoke was used to 
visualize the flow field in two ways. First, smoke was introduced in one tube 
to visualize the streak line. Next, a rake was used to visualize the entire 2-D 
flow field. 

Photos were taken using a Nikon 35mm camera and Kodak TMAX-400 
ASA black and white film. The film speed was set to 1/500 seconds with an 


aperture setting of 2.0 for the light conditions used at low rotational velocities. 
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At higher rotational velocities film speeds of 1/1000 and 1/2000 seconds were 


used with aperture settings of 2.0 and 1.4 respectively. 


E. RESULTS AND DISCUSSION 

The results for the first tunnel condition flow visualization experiment is 
shown in Figure 4.5. Here, the wind tunnel is turned off and the wave 
propeller is oscillating at approximately 1500 rpm. It was shown that the wave 
propeller accelerated the smoke through the tunnel and out the exit similar to 
a fan blade. This experiment was done for the cases where the wave propeller 
oscillated both clockwise and counterclockwise with the same result. This was 
clear evidence the wave propeller produced thrust; hence, showing that circular 
motion can simulate plunge motion which is known to produce thrust. 

The results of the second tunnel condition flow visualization experiment 
are shown in Figures 4.6 - 4.20. Figure 4.6 shows the wave propeller at 
steady state and O° AOA. It can be seen that the airfoil produces a highly 
turbulent wake at the very low Reynolds number (26,500) generated in this 
tunnel. The boundary layer would still remain attached, which was an 
improvement over the earlier airfoil design. Figures 4.7 through 4.14 show the 
streak line flow visualization experiment at various AOQA's, and frequencies. 
These photos show several important features about the street vortices 
produced by different oscillating conditions. Most of these pictures reveal the 


propulsive vortical street pattern discovered in reference 16. Unfortunately, the 
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wave propeller was constructed to oscillate at too high an amplitude (diameter) 
for this tunnel, and many of the vortical patterns would hit the walls and 
disperse. In Figure 4.10, the vortical pattern shows that the bottom vortex is 
rotating clockwise, and the top vortex is rotating counterclockwise, which is 
a thrust producing vortical street. It was observed that the AOA and frequency 
greatly affect the vortical strength (size). Increasing the AOA and frequency 
led to an increase in wake vorticity. Furthermore, these increases led to a 
stronger influence on the streak line in front of the wave propeller. Increasing 
the AOA beyond 5° led to stall, and subsequently, a decrease in thrust. Figure 
4.13 shows the wave propeller in a stalled condition (20° AOA), and the 
rotating vortices are both weaker and positioned such that they would not 
induce significant thrust. In fact, Figure 4.14 shows a dynamic stall vortex 
building up on the leading edge. 

Figure 4.15 - 4.18 contain the rake flow visualization experiment. Figure 
4.15 and 4.16 display the propulsive vortical pattern produced by the wave 
propeller for the 2-D flow field. Figures 4.17 and 4.18 indicate the large region 
of separation on the airfoil’s upper surface. Again, the dynamic stall vortex is 
seen to form at the leading edge in Figure 4.17. These pictures best reveal the 
stalled condition of the wave propeller operating at low Reynolds number and 
high amplitude. Both of these conditions were limitations in the design and test 


equipment. 
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Figure 4.5 Wave propeller operating at 1500 rpm, 0° AOA with the wind 
tunnel off. 








Figure 4.6 Steady state, O° AOA. 
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Figure 4.7 RPM = 1085, AOA = 0°, rotating counterclockwise. 
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Figure 4.8 RPM = 165, AOA = 5°, rotating counterclockwise. 


83 


ERR RRRVOL RR sce RR 





Figure 4.9 RPM = 620, AOA = 5’, rotating clockwise 
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Figure 4.10 RPM = 1085, AOA = 5°, rotating clockwise. 
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Figure 4.11 RPM = 620, AOA = 10°, rotating counterclockwise 





Figure 4.12 RPM = 1085, AOA = 10°, rotating counterclockwise 
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Figure 4.13. RPM = 620, AOA = 20’, rotating counterclockwise. 
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Figure 4.14 RPM = 1085, AOA = 20’, rotating counterclockwise. 
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Figure 4.15 RPM = 620, AOA = O°, rotating counterclockwise 





Figure 4.16 RPM = 1085, AOA = O°, rotating counterclockwise 








Figure 4.17 RPM = 1085, AOA = 0°, rotating counterclockwise 





Figure 4.18 RPM = 1085, AOA = O°, rotating counterclockwise 
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V. ENHANCED LIFT FLOW VISUALIZATION EXPERIMENT 


A. INTRODUCTION 

Experimental research in high angle-of-attack flight has long been an area 
of greatinterest. The area of active research centers around understanding and 
controlling the boundary layer to enhance lift in the post-stall region. To this 
end, many steady-state boundary layer control devices have been installed in 
an attempt to delay separation. 

Boundary layer control has also been attempted through unsteady 
excitation mechanisms. The most promising method Is the wave propeller first 
suggested by Wilhelm Schmidt in the German Journal of Flight Sciences in 
1965 [Ref. 14]. Schmidt’s wave propeller consisted of a single wing that 
performed a plunging motion perpendicular to the freestream. The oscillating 
wing was mounted between two stationary main lifting wings. In this 
arrangement the angle-of-attack of the wave propeller was varied to achieve 
optimum performance. Results from Schmidt's work indicate that with the 
optimum configuration the steady-state stall angle-of-attack could be increased 
to beyond 25° and the corresponding lift coefficient increased by four times. 

The results on flow visualization in this report was a continuation of the 
on-going research at the Naval Postgraduate School to better understand the 


beneficial effects of a wave propeller operating in the vicinity of a main wing. 
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The initial research and construction of the wave propeller was undertaken by 
a former NPS student, Carl Dane. Carl Dane’s thesis [Ref. 15] contains the 
results of pressure distributions on a main wing with the wave propeller 
configured to oscillate aft and slightly above the wing. Unfortunately, the wave 
propeller was designed poorly due to construction limitations and no beneficial 
effects were observed. The purpose of this research was to build a more 
efficient wave propeller; then, through flow visualization, understand the 


effects of the wave propeller on a stationary wing. 


B. EXPERIMENTAL SETUP 


1. Wave Propeller and Stationary Wing 

The wave propeller was the same as described in Chapter IV (Figure 
4.3). In this experiment, the wave propeller was mounted vertically in the NPS 
low speed wind tunnel behind a stationary wing. The wing hada 12-inch chord 
and a 28-inch span, the vertical length of the test section [Ref. 15]. The wing 
was constructed of a NACA 66(215)-216 airfoil section and could be set 
between O and 20 degrees AOA. Unfortunately, the wave propeller operated 
less efficiently in the vertical position. This was most likely a result of the 
increased weight on the bearings. A new power supply was used to provide 
increased rotational speeds of the wave propeller. A DC reversible power 
supply rated at a constant 48 volts, and 15 amps was used. With this power 


supply, the wave propeller operated at a rotational speed of 600 RPM. 
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2. Wind Tunnel 
The experiment was conducted in the 32x45 inch low speed wind 
tunnel located at the Navel Postgraduate School [Ref 17]. This tunnel was 
designed at the Aerolab Development Company of Pasadena, California as a 
closed circuit, single-return wind tunnel that is 64 feet long and 25.5 feet wide. 
The tunnel is powered by a 100 horse power motor which drives a three blade 
variable pitch fan by a four-speed International truck transmission. This motor 
and transmission can efficiently drive tunnel velocities up to 180 mph. Fan 
induced swirl is removed by eight stator blades located directly downstream of 
the fan. The tunnel cross section gradually expands from the fan to the settling 
chamber while turning through three sets of 90° corner vanes. The settling 
chamber includes two fine wire mesh anti-turbulence screens designed to break 
up large turbulent fluctuations. The flow accelerates to the test section though 
a 10:1 contraction cone. 
The test section is made up of transparent glass sidewalls (access doors) 
and upper wall for illumination and viewing. The test section is designed to 
operate at atmospheric pressure, and has a breather slot located aft of the test 


section to compensate for leakage losses throughout the tunnel. 


3. Smoke Generator 
The smoke was generated using the same Rosco smoke machine 
described in Chapter |V. The smoke was piped from the top of the tunnel into 


a blast tube. The blast tube directed the smoke into the tunnel slightly 
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upstream of the test section. In this configuration, the smoke was injected into 


the test section directly at the stationary wing. 


C. TEST PROCEDURE 

The wave propeller was mounted vertically from the top of the wind 
tunnel to extend down into the test section slightly behind the stationary wing 
(Figures 5.1 and 5.2). In this position, the wave propeller could rotate both 
clockwise or counterclockwise at various angle-of-attack settings. 

The tunnel was operated at speeds between 10-30 fps. The smoke was 
injected in the tunnel slightly upstream of the stationary wing. At this point, 
the wing angle-of-attack was increased until the boundary layer was just 
beginning to separate. Then, the wave propeller was turned on, and 
observations concerning the boundary layer were made. This continued for 


different AOA settings and rotational directions of the wave propeller. 


D. RESULTS AND DISCUSSION 

Figures 5.3 and 5.4 show the best pictures taken in the series of 
experiments. Figure 5.3 is a picture of the flow over the stationary wing with 
the wave propeller turned off. The flow is shown to be slightly separated. 
Figure 5.4 shows the flow over the stationary wing with the wave propeller on. 
Here, the flow is seen to be more attached. Although many pictures were 


taken, very few came out clearly because of the poor lighting conditions and 
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limited time available to set appropriate flow visualization conditions for this 
experiment. 

Visual observations of the wave propeller’s influence on the stationary 
wing were more convincing than the pictures. It was clear that the wave 
propeller had a beneficial effect on the stationary wing that energized the 
boundary layer at high angle-of-attack and delayed separation. It was observed 
that the rotation direction, AOA, and relative position of the wave propeller all 
controlled the degree of influence on the stationary wing. It was noted that a 
clockwise direction, and increase in AOA had a stronger positive influence on 
the stationary wing. It was unclear as to the best relative position of the wave 
propeller. This was complicated by the limited area available in the test section 
to move the wave propeller. It was further noted that the beneficial influence 
was Stronger at lower tunnel speeds. This was a result of the low rotational 
velocities achieved by the wave propeller. Apparently, the ratio of oscillation 
frequency to forward speed is critical for upstream influence. The higher the 
rotational speed of the wave propeller, the faster the freestream velocity can 


be and still achieve a beneficial upstream influence. 
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VI. CONCLUSION AND RECOMMENDATIONS 


A. SINGLE AIRFOIL ANALYSIS 

The non-linear theory presented here for harmonic motion, and the phase 
lag relationships that exist between the airfoil motion and build-up of 
aerodynamic forces has been extensively verified with existing linear theory. 
Furthermore, a non-linear theory has been presented that predicts the 
propulsive forces and efficiencies associated with pitching, and plunging 
airfoils. Finally, a theory has been presented that suggests a definite link 
between quasi-steady ramp motion (dynamic lift), and a delayed steady state 
event. 

More extensive studies of the quasi-steady phenomenon are recommended 
to thoroughly understand the relationships that exist between steady and 
unsteady motion. Another area of research is the study of the effect of airfoil 
geometry on aerodynamic flutter. A computational method can be developed 
using the presented codes and Theodorsen theory to solve the two-dimensional 


flutter determinant. 


B. TWO AIRFOIL ANALYSIS 
The modified version of USPOTF2A has been verified against existing 


theoretical studies. Furthermore, a limited theory has been presented that 
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predicts the thrust and efficiency associated with two airfoils in close 
proximity. Finally, it was shown in the upwind influence study that an 
oscillating airfoil will have a beneficial influence on the pressure distribution and 
pressure gradient of a stationary airfoil. 

In Chapter Ill the current limitations of the USPOTF2A code were noted. 
The limitation exists when the discrete vortices from one airfoil’s wake interact 
with the other airfoil. The code should be improved to eliminate this difficulty. 
This would greatly enhance the capability of this code when working with large 


amplitude motions or when the airfoils are in tandem. 


C. FLOW VISUALIZATION EXPERIMENTS 

The flow visualization experiment successfully showed the development 
of thrust produced by the wave propeller. Furthermore, the vortex street was 
analyzed to determine which conditions were more favorable for thrust 
production and which conditions induced stall. The enhanced lift flow 
visualization study was not a complete success. The wave propeller did not 
operate smoothly in the vertical position, and the upwind influence was only 
minor. Visual observation indicated a beneficial effect, but further photography 
is needed to document the event. Unfortunately, there was not sufficient time 
for experiments to produce adequate pictures of the phenomenon. 

It is recommended that further experiments be conducted in the low speed 


smoke tunnel. The wave propeller should be modified for a smaller rotational 
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diameter (smaller plunge amplitude). Smaller amplitudes will prevent vortex 
interactions with the tunnel walls. Finally, it is recommended that the wave 
propeller be rebuilt so binding in the bearings will not inhibit the motion in the 
vertical position. Further experiments with smoke visualization and pressure 
measurements are required to better understand the influence of the wave 
propeller on the stationary wing. Certainly, the ability to move the wave 
propeller to various positions around the stationary wing must be included in 


follow-on experiments. 
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APPENDIX A 


eCCC 
C THIS PROGRAM TAKES THE INPUT FILE (FILE CODE 14) CREATED BY 
C U2DIIF AND CONVERTS THE DATA TO A FREQUENCY, AMPLITUDE, AND 
C PHASE SHIFT. IT ALSO DETERMINES THE PROPULSIVE FORCES AND 
C EFFICIENCY. 
Eee 
PROGRAM PHASESHIFT 
DIMENSION PHASE(3) ,CL(400) ,CM(400) ,CD(400), 
+ ALPHA(400) , TIME(400) ,T(400) , AMP(3) ,CKT (400), 
+ FN(400) ,R(400) ,DAT(400) , FNT(3,400) , HY (400) 
mecucAl FLAG 
REAL L1,L2,L3,L4,M1,M2,M3,M4 
PI = ACOS(-1.0) 
C READ TYPE OF MOTION (0=PITCH, 1=PLUNGE) 
READ (14,*) MOTION, ALP1 
C READ NUMBER OF DATA POINTS AND FREQ 
READ (14,*) NPTS, W 
C READ DATA (TIME, CL,CM,CD) 
DO I=1,NPTS 
ReaD (s,*) TIMe (LI) ,CL(1L),CM(I),CD(T) 
END DO 
DO I=1,NPTS 
READ (12,*) T(1I),ALPHA(I) 
END DO 
DO I=1,NPTS 
READKT3,*) T(L) , HY (I) 
END DO 
me 200, J = 1,2 
DO I = 1,NPTS 
IF (J .EQ. 1) THEN 


DAT(I) = CL(T) 
ELSE IF(J .EQ. 2) THEN 
DAT(I) = CM(TI) 

END IF 


END DO 
C READ POSITION DATA 
IF (MOTION .EQ. 1) THEN 
DO I=1,NPTS 


ALPHA(I) = HY(T) 
END DO 
ZERO = .00001 
ELSE 
ZERO = .O1 
END IF 


ILS N=2 


SNS, 


DO f=2, NEES 
IF (ALPHA(I).LE.ZERO .AND. ALPHA(I) .GE. -ZERO) GOTO 16 
N=N+1 
END DO 
16 M=1 
DO I=N+1,NPTS 
IF (ALPHA(I) .LE.ZERO.AND.ALPHA(I) .GE. -ZERO) GO TO 10 
M = M+l1 
END DO 
ie. CALL AMPLITUDE (DAT, AMP, NPTS,J) 
PRINT*, ’AMPL = ’,AMP(J) 
@: 
C DETERMINE PHASE SHIFT 
é 
25 PHE==—0 
ERR .= T0000 
Oli = ppl 


PRINT *.> N. Mog S57 Neier 
ec 
C BEGIN ITTERATION TO CONVERGENCE 
& 
30 TE (ERR: {LbT.-TCOb) eco te.35 
SUM = 0 
DO I = N,N+M 
FN(I) = AMP(J)*SIN(W*T(I) + PHI) 
R(I) = ABS(FN(I) - DAT(I)) 
SUM = SUM + R(I) 
END DO 
IF((SUM - ERR) .GT. 0) THEN 
CN = CN/2.0 
NUM = NUM +1 
IF(NUM/2.0 .EQ. INT(NUM/2.0)) THEN 


FLAG = .FALSE. 
ELSE 
FLAG = .TRUE. 
END IF 
END IF 


IF (FLAG) THEN 
PHI = PHI + CN*PI/180.0 
COUNT = COUNT + 1 
ELSE 
PHI = PHI - CN*PI/180.0 
COUNT = COUNT + 1 
END IF 
IF (COUNT .GT. ITTER) GOTO .s00 
ERR = SUM 
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C 


35 


500 


100 


20 


200 


Gero 30 

PHASE(J) = PHI*180.0/PI 

Preeni*,’ ERR = ’, ERR 

Sono, 100 

PRINT*, ‘’ITTERATION GREATER THEN 500...INCREASING TOL’ 
TOL = 10*TOL 

COUNT = 0 

so, TO 30 

PF(J .EQ. 1) WRITE (6,20) PHASE(J)+180.0 

IF(J .EQ. 2) WRITE (6,21) PHASE(J) 


FORMAT(1X,’ CL PHASE =',F12.6) 
FORMAT (1X,' CM PHASE =',F12.6) 
DO I = 4,NPTS 
FNT(J,1I) = AMP(J)*SIN(W*T(I) + PHASE(J) *PI/180.0) 
END DO 
CONTINUE 


DO I = 4,NPTS 
Miembro j= Je (h) PNT (la); PNG (2,1) jee NT (3,2) 
END DO 


C DETERMINE THE PROPULSIVE EFFICIENCY 


c 


ee” 
Il 


ay PHASE (1) *PI/180.0 
PHASE(2) = PHASE(2)*PI/180.0 


DO I =2,N+M 
Cblor = CDTOT@ CD(1I) = €d(1) 
K = K+l 
END DO 
DBAR = CDTOT/K 
DBAR = DBAR 
PRINT*, ’AVERAGE DRAG (THRUST) , TOTAL DRAG (THRUST) 
OB AR CDOTOT 
IF(MOTION .EQ. 1) THEN 


WBAR = -.5*W*SIN(PHASE(1) ) *AMP (1) 
ETAS = 2*DBAR/WBAR 

ELSE 
WBAR = .5*W*SIN(PHASE (2) ) *AMP(2) 
ETAS = DBAR/WBAR 

END IF 

PRINT*,’ETAS,WBAR = ’,ETAS,WBAR 


SUBROUTINE AMPLITUDE (DAT, AMP,NPTS,J) 
DIMENSION DAT (200) ,AMP(3) ,AMP1(10) , AMP2 (10) 


N2 = 0 

M2 = 0 

WO T= 1,10 
AMP1(I) = 0 
AMP2 (I) = 0 
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END DO 
DO £ = 2 ,NPrs-a 
IF(DAT(I+1) .LT. 0 .AND. DAT(L) (aitGw eon 
IF (ABS (DAT(I+1)) .GT.ABS(DAT(I))) THEN 
IF((N2+1)/2.0 .EQ. INT((N2+1)/2.0)) N2 = sie 
TMP = ABS (DAT (I+1) ) 
ELSE 
GOTO 10 
END IF 
IF(TMP .GT. AMP1(N2)) THEN 
AMP1(N2) = TMP 
TMP =0 
END IF 
ELSE 
TF ((N2+1)/2.0 .NE. INT((N2+1)/2.0)) WN2 Sone 
IF (ABS (DAT(I+1)) .GT.ABS(DAT(I))) THEN 
TMP = ABS (DAT(I+1) ) 
ELSE 
GOnos 10 
END IF 
IF(TMP .GT. AMP2(N2)) THEN 
AMP2(N2) = TMP 
TMP =0 
END IF 
END IF 
10 END DO 
IF (AMP1(2) .GT. AMP2(2)) THEN 
IF(AMP1(2) .LT. AMP2(3)) THEN 
COMP = AMP1(2) 
ELSE 
COMP = AMP2(3) 
END IF 
ELSE 
IF(AMP2(2) .LT. AMP1(3)) THEN 
COMP = AMP2 (2) 


ELSE 
COMP = AMP1 (3) 
END IF 
END IF 
PRINT*,’COMP = ’,COMP 


DO eile 224 NZ 
IF (ABS (AMP1(1I)-COMP) .GT. .1*COMP) GO TO 20 
M2 = M2 +1 


AMP(J) = AMP(J) +AMP1(T) 
GO TO. 30 
20 IF (ABS (AMP2 (I) -COMP) .GT. .1*COMP) GO TO 30 
M2 = M2 +1 
AMP (J) = AMP(J) + AMP2(T) 
30) PRINT*,’AMP(I), AMP(I) = ’, AMP1(I), AMP2(I) ,M2 
END DO 
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IF (J .LT. 3) THEN 
AMP (J) = AMP(J) / (M2) 
ELSE 
AMP (J) = AMP(J) /(2*M2) 
END IF 
RETURN 
END 
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APPENDIX B 


CeCe 
C THIS PROGRAM TAKES THE INPUT FILE (FILE CODE 14) CREATED BY 
C USPOTF2A AND CONVERTS THE DATA TO A FREQUENCY, AMPLITUDE, 
C AND PHASE SHIFT. IT ALSO COMPUTES THE PROPULSIVE FORCES AND 
C EFFICIENCIES ASSOCIATED WITH THAT DATA SET. 
CCCE 
PROGRAM PHASESHIFT2 
DIMENSION PHASE (2,2) ,CL(2,450) ,CM(2, 450) , CD(2) 4508 
+ K(2) ,ALPHA(2,450) ,TIME(450) ,T(450) ,AMP(2,2) ,ALP(2), 
+ CDTOT(2),FN(400) ,R(450) ,W(2) ,DAT(2,450) , FNT(2,450,2), 
+ HY(2,450) ,T1 (450) ,CP1(450) ,CP2 (450) ,CP3 (450) , CE4aw 
REAL L1,L2,L3,L4,M1,M2,M3 ,M4 
LOGICAL FLAG, FIRST 
Pl = ACOS(-1-. 0) 
C READ TYPE OF MOTION (0=PITCH, 1=PLUNGE) 
READ (14,*) MOTION, ALP1 
C READ NUMBER OF DATA POINTS AND FREQ 
READ(14,*) NPTS, W(1),W(2) 
PRINT*, MOTION, ALP1 
PRENE* ,NE@S, Wi lic, WZ) 
GC READ DATA. (EIME, Ch,CM Cr) 
READ (7, *) 
READ (7, *) 
READ (7, *) 
DO I=1,NPTS 
READ (7,*) TIME(I),CL(1,I),CL(2,1),CM(1, 1) Gia 
+ Coy CDi 2 4) 
END DO 
DO I=1,NPTS 
READ (12,*) T(I) ,ALPHA(1,I) ,ALPHA(2,T) 
END DO 
PO f=1, Nels 
READ(13,*) T(1), AY (31) ree 


END DO 
IF (MOTION .EQ. 1) THEN 
VO? Lig=mre 2 
DO I=1,NPTS 
ALPHA(L,I) = HY(L,I) 
END DO 
END DO 
ZERO = .00001 
ELSE 
ZERO = .01 
END IF 


IF (W(1) .EQ. 0) THEN 
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FIRST = .FALSE. 
Wwe = W(2) 
DO I = 1,NPTS 
ALPHA(1,I) = ALPHA(2,TI) 
END DO 
ELSE IF (W(2) .EQ. 0) THEN 
FIRST =.TRUE. 
W(2) = W(1) 
DO I = 1, NPTS 
ALPHA(2,I) = ALPHA(1,1TI) 
END DO 
END IF 
Be 200) ime 1, 2 
Oot — |, NPTS 
IF (J .EQ. 1) THEN 
Air) ae= Cla aL) 
PAM G2) Pie -vCh(2.. 8) 
ELSE IF(J .EQ. 2) THEN 
BAT (1, F) 5= CM(1, 1) 
Pek ee CM (2 Dh) 
END IF 
END DO 
C READ POSITION DATA 
BOsrLOG tf = 1,2 
Ls N=2 
DO I=2,NPTS 
IF(ALPHA(L,I) .LE. ZERO .AND. ALPHA(L,I) .GE. 
+ -~ZERO) GO TO 16 
N=N+1 
END DO 
16 M=1 
DO I=N+1,NPTS 
IF (ALPHA(L,I) .LE. ZERO .AND. ALPHA(L,I) .GE. - ZERO) 
+ GO TO 10 
M = M+1 
END DO 
10 PRINT*,W(L) 
CALL AMPLITUDE (DAT, AMP,NPTS,L,J,N) 
PRINT*, ‘’AMPL = ’,AMP(J,L) 
& 


C DETERMINE PHASE SHIFT 
c 


25 PHI = 
ERR = 10000 
TO -= .01 
CN = 1.0 
ITTER = 500 
COUNT = 0 
NUM = 1 
FLAG = .TRUE 
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PRINT* ,9N,M,b,0 =) 2 Nye 


C BEGIN ITTERATION TO CONVERGENCE 


30 


Sys) 


SOC 


100 


20 


ext 


IF({ERR .LT. TOL) Gonteess 


SUM = 0 

DO I = N,N+M 
FN(I) = AMP(J,L) *SIN(W(L) *T(I) + PHI) 
R(I) = ABS(FN(I) - DAT(L,I)) 
SUM = SUM + R(I) 

END DO 


IF( (SUM — ERR) {Gi 0 0) aera 
CN = CN/7Z0 
NUM = NUM + 1 
IF (NUM/2.0 .EQ. INT(NUM/2.0)) THEN 


FLAG = .FALSE. 
ELSE 

FLAG = .TRUE. 
END IF 


END IF 
IF (FLAG) THEN 
PHI = PHI + CN*PI/180.0 
COUNT = COUNT + 1 
ELSE 
PHI = PHI - CN*PI/180.0 
COUNT = COUNT + 1 

END IF 

IF (COUNT .GT. ITTER) GOTO 500 

ERR = SUM 

GOTO 30 

PHASE(J,L) = PHI*180/PI 

PRINT*, ‘PHASE(J,L) = ’, PHASE(d, Lb) 

PRIENT*. “BRROR =~. VERE 

GO TO 100 

PRINT*, ’ITTERATION GREATER THEN 500...INCREASING TOL’ 
TOL = 10*70b 

CN = 10*CN 

COUNT = O 

GO TO 30 
CONTINUE 

IF(J .EQ. 1) WRITE (6,20) PHASE(J,1)+180.0 ,PHASE(J,2) 
+180.0 
IF(J .EQ. 2) WRITE (6,21) PHASE(J,1), PHASE (Ge 


FORMAT (1X,’ CL PHASE(1) =’ ,F10.6/)52%, “ CONPHASr 
=? 1) 202.6) 
FORMAT (1X,' CM PHASE(1) =’,F10.6,5X, CM PHASa 
=" 30. 6) 


DO I = 4,NPTS 

FNT (J,1I,1)=AMP(J,1) *SIN(W(1) *T(I)+ PHASE(J,1)*PI/180.0) 
FNT(J,1I,2) =AMP(J,2) *SIN(W(2) *T(I)+ PHASE(J,2) *PI/180.0) 
END DO 

PHASE (J,1) = PHASE(J,1) *P1/180.0 
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PHASE(J,2) = PHASE(J,2) *PI/180.0 
200 CONTINUE 
c 
C DETERMINE THE PROPULSIVE EFFICIENCY 
6 
SeLNPUT DTS 
PRINT*, ’INPUT DTS ' 


READ*, DTS 
Kl = Q 
(Orin — ay, 2 
CDTOT(L) = 0 
re) = 0 
Ow =DTS, 2*DTS 
CDTOT(L) = CDTOT(L) + CD(L,TI) 
K(L) = K(L) + 1 
END DO 
END DO 
DBAR1 = CDTOT(1) /K(1) 
DBAR2 = CDTOT(2) /K(2) 


PRINT*,’AVERAGE DRAG, TOTAL DRAG (AF 1) 
ae DBARI CDTOT (1) ,K(1) 
PRINT*, ’AVERAGE DRAG, TOTAL DRAG (AF 2) 
Pee BAKZ, COTOT(2) , K(2) 
Pr (MOTION .EO. 1) THEN 
PieiRST) THEN 


PHAS = PHASE(1,1) 
AMPL = AMP(1,1) 
ELSE 
PHAS = PHASE(1, 2) 
AMPL = AMP (1,2) 
END IF 
WBAR = -.5*W(1)*SIN(PHAS) *AMPL 
ELSE 
iP (FIRST) THEN 
PHAS = PHASE (2,1) 
AMPL = AMP(2,1) 
ELSE 
PHAS = PHASE (2,2) 
AMPL = AMP (2,2) 
END IF 
ALP1 = ALP1*PI/180.0 
WBAR = -.5*W(1)*SIN(PHAS) *AMPL 
END IF 
ETAS = 2* (DBAR1+DBAR2) /WBAR 
PRINT*, 'ETAS,WBAR = ',ETAS,WBAR 


DO I=DTS,2*DTS 
Chinon = ChLirorT + CL(1,1) 


CEZLor CLZrToO Ta Cli2 1) 
Role = a) 
END DO 
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CLIAVG = Chileon ta 

CL2AVG = CL2TOT/K1 

DELCL1 = CLIAVG - CL(1,1) 

DELCL2Z = CL2ZAVG seni.) 

PRINT*, ’CLIAVG, DELCL1 = “ ,ChLIAVG PE iere 
PRINT*, ’CL2AVG, DELCL2 = ’ , Ch2AVG 7 PEEeEZ 


DO I = 4,NPTS 
WRITE (15,*)T(I),FNT(1,1,1),FNT(1,1,2), PNI(2 eee 
+ ENE (2; L332) 
END DO 


END 


SUBROUTINE AMPLITUDE (DAT, AMP,NPTS,L,J,N) 
DIMENSION DAT(2,200) ,AMP(2,2),AMP1(10) , AMP2 (10) 


N2 = 0 

M2 = 0 

DOF i=) 20 
AMP1(I) = 0 
AMP2 (I) = 0 

END DO 


BPOsE = 2, NETS-1 
IF (DAT(L,I+1) .LT. 0 .AND. DAT(L,1I) .LT. 0) ) etna 


IF (ABS (DAT(L, 1I+1)) .GT. ABS(DAT(L,1)}) THEM 
IF((N2+1)/2.0 .EQ. INT((N24+1)/2.0)) N2 = Noa 
TMP = ABS (DAT(L,1+1) ) 
ELSE 
GOTO 10 
END IF 


IF(TMP .GT. AMP1(N2)) THEN 
AMP1(N2) = TMP 
TMP =0 
END IF 
ELSE 
IF((N2+1)/2.0 .NE. INT((N2+1)/2.0)) N2 =e 
IF (ABS (DAT(L,I+1)) .GT.ABS(DAT(L,I))) THEN 
TMP = ABS (DAT(L,I+1) ) 
ELSE 
GOTO 10 
END IF 
IF(TMP .GT. AMP2(N2)) THEN 
AMP2 (N2) = TMP 
TMP =0 
END IF 
END IF 
10 END DO 
IF (AMP1(2) .GT. AMP2(2)) THEN 
IF(AMP1(2) .LT. AMP2(3)) THEN 
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12, 


20 


al. 


30 


COMP = AMP1(2) 
ELSE 

COMP = AMP2 (3) 
END IF 


ELSE 
IF(AMP2(2) .LT. AMP1(3)) THEN 
COMP = AMP2(2) 


ELSE 
COMP = AMP1 (3) 
END IF 
END IF 
PRINT*,’COMP = ’,COMP 


DO Gie= 2,N2 
IF (AMP1(I) .GT. COMP) GO TO 19 
IF (ABS (AMP1(I)-COMP) .GT. .12*COMP) GO TO 20 
M2 = M2 +1 
AMP(J,L) = AMP(J,L) +AMP1 (I) 
GO TO 30 
IF (AMP2 (I) .GT. COMP) GO TO 21 
TF (ABS (AMP2(I)-COMP) .GT. .12*COMP) GO TO 30 


M2 = M2 +1 
AMP(J,L) = AMP(J,L) + AMP2(T) 
PRINT*, ’AMP(I), AMP(I) = ’, AMP1(I), AMP2(I) ,M2 
END DO 
AMP(J,L) = AMP(J,L) / (M2) 
RETURN 


END 
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